#!/usr/bin/Rscript

#plotEStationData.R:
#
# Last modified 14 August 2008
# Author: Cameron Bracken
#         Student Researcher
#         NOAA National Weather Service
#         300 Startare Drive
#         Eureka, CA 95501
#         cameron.bracken@noaa.gov
#
#
#

fort.61 = '../adcirc/fort.61'
outfile = 'plots/elevation.pdf'


    #read adcirc data 
con = file(fort.61)

stat = scan(fort.61,skip=1,nlines=1,quiet=T)
ntimes = stat[1]
nstations = stat[2]
adcElev = matrix(0,nrow=ntimes,ncol=nstations)

adctimes = 0
skip = 2 

for(i in 1:ntimes){
    adctimes[i] = scan(con,skip=skip,nlines = 1,quiet=T)[1]
    adcElev[i,] = scan(con,skip=(skip+1),nlines=nstations,quiet=T)[seq(2,2*nstations,by=2)]
    skip = skip + 1 + nstations
}
close(con) 
adctimes=as.POSIXct(adctimes,tz="UTC",origin="2004-01-01")


stimestamp = format(range(adctimes)[1],'%Y%m%d')
etimestamp = format(range(adctimes)[2],'%Y%m%d')

pdf(outfile,width = 18, height = 4)  


    #download and parse co-ops elevation data 
    
        #make sure times are in local time
    adctimes = c(adctimes)


        #Now plot it
for(i in 1:nstations){
    ylim = c(min(c(adcElev[,i])),max(c(adcElev[,i])))
    plot(adctimes,adcElev[,i],col='dodgerblue',type='l',ylab='Water Level',xaxt='n')
    #legend('topright',c('Observed','Adcirc Predicted'),col=c(1,'dodgerblue'),lty='solid')

    r <- as.POSIXct(round(range(c(adctimes)), "days"),tz="PST")
    axis.POSIXct(1,r, at=seq(r[1], r[2], by="days"), format="%b %d")
}
dev.off()

st <- as.integer(substr(stimestamp,1,4))
st[2] <- as.integer(format(adctimes[1],"%j"))*24*4
print(st)
f <- 35040
adcirc.ts <- ts(adcElev[,1],start=st,frequency=f)

write(t(cbind( format(adctimes,"%Y,%m,%d,%H,%M"),adcElev[,1])),sep=',',file='adcirc-ts.out',ncolumns=2)
save(adcirc.ts,file='adcirc-ts.RData')
cat('Wrote elevation time series to adcirc-ts.out and adcirc-ts.RData\n')
